      subroutine pwrspec(nx,dt,data,freq,pwr)
      dimension data(nx),freq(nx),pwr(nx)
 
      kmax = 20
      nexp=0
      ntmp=nx
 
      do j = 1,kmax
      ntmp=ntmp/2
      if(ntmp.eq.1)go to 101
      nexp=nexp+1
      end do
101   continue
      npad =  2**(nexp+1) - nx 
c     print *,'npad,nx ',npad,nx
c     print *,'2^n ',2**nexp
      


      if(npad.ne.0)then
c      
c      pad data with zeros to power of 2
c
       do 21 j=nx+1,2**(nexp+1)
       data(j)=0.
21     continue
       nx=2**(nexp+1)
       else
       endif

c      print *,'nexp ',nexp
c      do j=1,nx
c      print *,data(j)
c      end do
 
       nby2=nx/2

      do j = 1,nby2
      freq(j)=0.
      pwr(j)=0.
      end do 
 
      fc=1./(2.*dt)
      isign=1
      call realft(data,nby2,isign)          
	freq(1)=0.
      pwr(1)=data(1)**2/(nx*nx)
      do 11 k=2,nby2-1
	freq(k)=float(k-1)*fc/float(nx)
	pwr(k)=(data(k+1)**2+data(nx-(k+1))**2)/(nx*nx)
11    continue
	freq(nby2)=.5*fc
      pwr(nby2)=data(2)**2/(nx*nx)
      return
      end

      SUBROUTINE REALFT(DATA,N,ISIGN)
      DIMENSION DATA(*)
      THETA=6.283185307/(2.0*N)
      C1=0.5
      IF (ISIGN.EQ.1) THEN
        C2=-0.5
        CALL FOUR1(DATA,N,+1)
      ELSE
        C2=0.5
        THETA=-THETA
      ENDIF
      WPR=-2.0*SIN(0.5*THETA)**2
      WPI=SIN(THETA)
      WR=1.0+WPR
      WI=WPI
      N2P3=2*N+3
      DO 11 I=2,N/2+1
        I1=2*I-1
        I2=I1+1
        I3=N2P3-I2
        I4=I3+1
        WRS=(WR)
        WIS=(WI)
        H1R=C1*(DATA(I1)+DATA(I3))
        H1I=C1*(DATA(I2)-DATA(I4))
        H2R=-C2*(DATA(I2)+DATA(I4))
        H2I=C2*(DATA(I1)-DATA(I3))
        DATA(I1)=H1R+WRS*H2R-WIS*H2I
        DATA(I2)=H1I+WRS*H2I+WIS*H2R
        DATA(I3)=H1R-WRS*H2R+WIS*H2I
        DATA(I4)=-H1I+WRS*H2I+WIS*H2R
        WTEMP=WR
        WR=WR*WPR-WI*WPI+WR
        WI=WI*WPR+WTEMP*WPI+WI
11    CONTINUE
      IF (ISIGN.EQ.1) THEN
        H1R=DATA(1)
        DATA(1)=H1R+DATA(2)
        DATA(2)=H1R-DATA(2)
      ELSE
        H1R=DATA(1)
        DATA(1)=C1*(H1R+DATA(2))
        DATA(2)=C1*(H1R-DATA(2))
        CALL FOUR1(DATA,N,-1)
      ENDIF
      RETURN
      END
      SUBROUTINE FOUR1(DATA,NN,ISIGN)
      DIMENSION DATA(*)
      N=2*NN
      J=1
      DO 11 I=1,N,2
        IF(J.GT.I)THEN
          TEMPR=DATA(J)
          TEMPI=DATA(J+1)
          DATA(J)=DATA(I)
          DATA(J+1)=DATA(I+1)
          DATA(I)=TEMPR
          DATA(I+1)=TEMPI
        ENDIF
        M=N/2
1       IF ((M.GE.2).AND.(J.GT.M)) THEN
          J=J-M
          M=M/2
        GO TO 1
        ENDIF
        J=J+M
11    CONTINUE
      MMAX=2
2     IF (N.GT.MMAX) THEN
        ISTEP=2*MMAX
        THETA=6.2831853/(ISIGN*MMAX)
        WPR=-2.0*SIN(0.5*THETA)**2
        WPI=SIN(THETA)
        WR=1.0
        WI=0.0
        DO 13 M=1,MMAX,2
          DO 12 I=M,N,ISTEP
            J=I+MMAX
            TEMPR=(WR)*DATA(J)-(WI)*DATA(J+1)
            TEMPI=(WR)*DATA(J+1)+(WI)*DATA(J)
            DATA(J)=DATA(I)-TEMPR
            DATA(J+1)=DATA(I+1)-TEMPI
            DATA(I)=DATA(I)+TEMPR
            DATA(I+1)=DATA(I+1)+TEMPI
12        CONTINUE
          WTEMP=WR
          WR=WR*WPR-WI*WPI+WR
          WI=WI*WPR+WTEMP*WPI+WI
13      CONTINUE
        MMAX=ISTEP
      GO TO 2
      ENDIF
      RETURN
      END
